A Monte-Carlo model for elastic energy loss in a hydrodynamical background 
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We present a computation of elastic energy loss of hard partons traversing the bulk hydrodynam- 
ical medium created in ultrarelativistic heavy-ion collisions. The model is based on perturbative 
Quantum Chromodynamics (pQCD) cross sections for 2 — > 2 processes in which a hard incoming 
parton is assumed to interact with a thermal parton from the medium. We model the interactions 
of this type in a Monte-Carlo framework to account properly for exact energy-momentum conser- 
vation, non-eikonal parton propagation, parton conversion reactions and the possibility to create 
additional hard recoiling partons from the medium. For the thermodynamical properties of the 
medium we use a hydrodynamical evolution model. We do not aim at a full description of high 
transverse momentum (Pr) observables at this point. Rather, we view the model as a starting point 
in obtaining a baseline of what to expect under the assumptions that the medium is describable 
by thermal quasifree partons and that their pQCD interactions with the high-energy partons are 
independent. Deviations from this baseline then call for more sophisticated medium description, as 
well as inclusion of higher-order processes and coherence effects in the pQCD scatterings. 

PACS numbers: 25.75.-q, 25.75.Bh 



I. INTRODUCTION 

The energy loss of hard partons propagating through the 
soft medium created in heavy-ion collisions has long been 
regarded as a promising tool to gain information on the 
medium properties [lMg]. Radiative energy loss, i.e. the 
idea that medium-induced radiation predominantly car- 
ries away energy from a hard parent parton has been 
rather successful in describing not only the hadronic nu- 
clear suppression factor for central collisions but also the 
effects of changing medium geometry Q. Calculations 
within dynamical evolution models in various formalisms 
fcl-fl0| have improved on this result and show also agree- 
ment with measured hard back-to-back correlations [TTI - 
Ojl and the measured suppression of protons 

Nevertheless, there are indications that a radiative en- 
ergy loss picture fails to describe the suppression of heavy 
quarks as seen in the measurements of the single-electron 
spectra llal. Energy loss due to elastic collisions with the 
medium [16-24J has been suggested as a possible solution 
to this problem. Such calculations indicate a large com- 
ponent of elastic energy loss also for light quarks and 
gluons. However, the pathlength dependence of elastic 
energy loss does not agree with the observed response of 
the suppression to a large pathlength bias as probed in 
back-to-back hadron correlations |25j . 

In the light of these inconclusive findings, it would be de- 
sirable to have a baseline calculation of elastic energy loss 
which treats the relevant physics in as detailed manner 
as possible. In particular, many analytical calculations of 
elastic energy loss assume eikonal propagation of a hard 
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parton, however in a more realistic description even hard 
partons, pr ^ 5 GeV, can be deflected from their path 
by interactions with the medium. Similarly, in processes 
like gg — > qq there is not only a change of energy between 
incoming and outgoing parton but also a change in par- 
ton species which influences the subsequent interaction 
probability of the outgoing parton. Finally, in scatter- 
ing processes in which the energy is distributed almost 
evenly among the outgoing partons, it is not justified to 
use the picture of energy from a hard parton being lost 
into the medium. Clearly, one also has to account for sit- 
uations in which both partons after a scattering process 
are hard or both become part of the medium, or in which 
the identity of the leading parton changes. 

In the present manuscript, we present work towards such 
a baseline computation. For this purpose, we develop 
a Monte Carlo (MC) simulation for the hard parton's 
interaction with the medium, as this is an appropriate 
framework to account for all the above mentioned ef- 
fects. In addition, we also try to to give estimates on 
the uncertainties of such a computation. Similar, in 
some aspects even more ambitious approaches are the 
MC models JEWEL (Jet Evolution With Energy Loss) 
26] and MARTINI (Modular Algorithm for Relativistic 
Treatment of heavy IoN Interactions) 27], both of which 
consider pQCD interactions of a parton shower with a 
thermal medium and include both elastic and radiative 
energy-loss components. In JEWEL, one models shower 
evolution and medium-induced radiation in a static (no 
flow, T constant) medium, and in MARTINI the vacuum 
showering is assumed to cease when the hydrodynamic 
evolution starts, while our approach focuses on the lead- 
ing partons. Like MARTINI, also we utilize a realistic 
well-tested 3-dimensional hydrodynamical description of 
the medium (in our case azimuthal symmetry in central 
collisions, and longitudinal boost symmetry but 3-d flow) 
instead of a static medium applied in JEWEL. 
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II. THE MODEL 

In the following, we make the assumption that the 
medium can be characterized as a thermal gas of quasi- 
free quarks and gluons and that its interactions with a 
hard parton can be treated as incoherent 2 — > 2 processes 
in leading-order pQCD. We stress that these assumptions 
need not be realized for the medium produced in ultrarel- 
ativistic heavy-ion collisions — any discrepancy between 
the calculation presented here and the data will reflect 
this. 

We model the medium by means of a hydrodynamical 
calculation [28| with the initial state obtained from the 
saturation of the phase space with minijets 29]. This 
provides the medium temperature T at each spacetime 
point. For the moment, we assume that only the scat- 
terings with a partonic medium are relevant and ignore 
the possibility that a hadronic medium might interact 
with a hard parton. Thus the mean free path of hard 
partons becomes infinite once the medium temperature 
drops below the critical temperature Tq- Given the lo- 
cal T, the distribution of scattering partners is given by 
the Fermi-Dirac (Bose-Einstein) distribution for quarks 
(gluons) respectively. We also account for the 3-d hydro- 
dynamical flow of the fluid cell. 

Our approach, most specifically the assumption that the 
relevant interactions can be described in terms of par- 
tonic pQCD interaction, shares many features with par- 
ton cascade models (e.g. (30l - [32l |). However, unlike in 
parton cascade models where the partonic pQCD inter- 
actions are assumed to describe the whole off-equilibrium 
bulk medium, we assume that the bulk is essentially ther- 
malized and use explicit pQCD interactions only to ac- 
count for the non-equilibrium dynamics between hard 
parton and medium. 

In the following, we first give a more detailed account of 
the hydrodynamical simulation, the parton-medium scat- 
tering rates and details of the MC implementation. Then 
we present our results for the partonic and hadronic nu- 
clear modification factor Raa{pt) and for QCD-matter 
tomography, and also discuss the relevant uncertainties. 
Some technical details regarding the scattering rate cal- 
culations can be found in the Appendix. 

A. Hydrodynamics 

In the present exploratory stage we consider only the 
cleanest possible case, (nearly) central A + A collisions 
in the mid-rapidity region. The produced QCD matter, 
which evolves and flows 3-dimensionally, can to a first 
realistic approximation then be described by using the 
nonviscous, azimuthally symmetric, longitudinally boost- 
invariant hydrodynamic framework. The 1+1 d hydro- 
dynamical framework, which we apply here for the QCD 
medium evolution, is presented in detail in Ref. [28| . 
The initial energy density e(r, tq) — eo(r) for the hy- 
drodynamics is here obtained from pQCD-minijet pro- 



duction [33| and final-state saturation as computed in 
the EKRT-model (34[. The transverse profile of eo is 
assumed to scale with the number of binary collisions, 
and thermalization is assumed to take place at forma- 
tion (i.e. at saturation) . For central Au-I- Au collisions at 
V$nn ~ 200 GeV, we thus initialize the hydrodynamics 
at (see Table 1 in [28]) r — tq = 0.17 fm/c, where the 
longitudinal proper time is defined as r = \Jt 2 — z 1 . The 
transverse flow is assumed to be zero at To. The Equa- 
tion of State (EoS) in our hydrodynamic background 
contains a high-T QGP phase, a low-T hadron reso- 
nance gas (HRG) phase, and a mixed phase between 
them. The QGP here is an ideal gas of gluons and 
three flavours of massless quarks and antiquarks obeying 
a Bag-EoS, while the HRG is an ideal gas of all hadrons 
and hadron resonances up to M — 2 GeV. The Maxwell 
construction gives the mixed phase, and the Bag constant 
B 1 / 4 = 239 MeV fixes the phase transition temperature 
to T c = 165 MeV. We terminate the hydrodynamic evolu- 
tion at a single decoupling temperature Td ec = 150 MeV. 
As discussed in Ref. [28j . the measured pion and kaon 
spectra and proton multiplicity in central Au+Au colli- 
sions at RHIC can be reproduced very nicely with such 
hydrodynamical set-up, once the thermal hadron spectra 
have been computed via the Cooper- Frye decoupling pro- 
cedure and once the strong and electromagnetic hadronic 
decays have been accounted for. For more details, as well 
as for plots of the initial energy density and of the phase 
boundary locations in the (r, r) plane, we refer the reader 
to Ref. g|. 



B. Scattering rate and mean free path 

Our simulation of energy losses of high-energy partons in 
the produced QCD matter rests largely on one particular 
physical quantity, namely the scattering rate for a high- 
energy parton of a type i, 

T l (p 1 ,u(x),T(x)) = ^2 r iJ ^ k i(p 1 ,u(x),T(x)), (1) 

where all possible elastic and inelastic partonic processes 
ij — > kl are accounted for by summing over all types 
of collision partners, j = it, d, s, u, d, s, g in the initial 
state, and over all possible parton type pairs (kl) in the 
final state. In general, the scattering rate depends on 
the frame, and in particular on the high-energy parton's 
4-momcntum pi, on the flow 4- velocity u(x) and on the 
temperature T(x) of the fluid at each space-time location 
x. In the notation below, we leave these dependences 
implicit but in the simulation they are fully accounted 
for. Ignoring Pauli blocking and Bose enhancement, we 
express the scattering rate as follows (see also Ref. |35j ]) 
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state particles are accounted for, the cross section is 



d 3 p 2 



d 3 P3 



2E 1 J (2tt) 3 2E 2 J (2n) 3 2E 3 J (2T:) 3 2E i 
fj(P2 ■ u,T) \M\^ k i(s,t,u) S 2 (s,t,u) 
(2n) i S^(p 1 +p 2 -p 3 -p 4 ). 



(2) 



where \M\f^ kl is the spin- and colour-summed/averaged 
scattering amplitude which depends on the standard 
Mandelstam variables s, t and u. The distribution func- 
tion for a thermal particle of momentum p 2 is denoted 
by fj(p 2 • u,T), which is the Bose-Einstein distribution 
f g = /s(p2 • U ,T) for gluons and the Fermi-Dirac dis- 
tribution f q = /d(p2 ■ u,T) for quarks. Since the initial 
state spins and colours have been averaged over in the 
squared matrix elements, the spin- and color-degeneracy 
factors in Jb and fjj are the standard ones, g g = 2-8 = 16 
for gluons, and g q — 2 ■ 3 = 6 for quarks and antiquarks, 
while the possible flavour-degeneracy needs to be deter- 
mined separately in each process. In our simulation, we 
evaluate the scattering rate always in the local rest-frame 
of the fluid, so that p 2 ■ u = E 2 and E\ is the energy of 
the high-energy parton i in this frame. 

The singularities arising in forward and backward scat- 
terings of massless partons — appearing in the scatter- 
ing amplitudes at t, u — > — need to be regularized 
in some physically meaningful way. For the moment, 
an exact field-theoretical solution to this non-equilibrium 
problem is not known. One possible treatment is to ap- 
ply resummed thermal (equilibrium) propagators sepa- 
rately for the longitudinal and transverse components, 
as is done in the AMY approach [U, |3(| . Since the prob- 
lem is, however, still phenomenological, we choose here 
a more transparent approach (tested also in JEWEL be- 
fore) and introduce a thermal-mass-like overall cut-off 
scale to = s m g s T for the momentum phase space, such 
that u,t < —to 2 . Here g s is the strong coupling constant 
and s m is a parameter of the order of one but which we 

2 

can vary. As we cannot control the running of a s = f*- , 
we keep the strong coupling constant fixed with momen- 
tum scale but vary both g s and s m in different runs of 
the MC code to test the model sensitivity to these pa- 
rameters. Since s + u + t = 0, the cut-off conditions are 
equivalent to — s + to 2 < t < — m 2 with the requirement 
s > 2m 2 . These are encoded in the function S 2 above, 
which is defined as 

S 2 (s, t, u) = 6(s > 2m 2 )6(-s + m 2 <t< -to 2 ) (3) 



For the evaluation of the Lorcntz-invariant part in 
Eq. (J2j) , it is useful to recall that the scattering rate can 
be expressed in terms of the scattering cross sections. For 
the process ij — > kl, where the Lorentz-invariant cut-off 
conditions are imposed, and the identities of the final 



ij—thl 



00 =- 



1 1 



d 3 P3 



d 3 p 4 



-S 2 (s,t,u) 



1 + S k i 2s J (2tt) 3 2E 3 J (2tt) 3 2E. 
■ (2tt)H (4 Hpi +P2-P3-P4) \M\ 2 3 _ kl (s,t,u), 

2 

1 1 



1 + 5 k l 167TS 2 



s+ra 2 



dt\M\%^ M , 



and we can rewrite Eq. $2$ as 
d 3 p 2 



ij—¥kl 



f j (E 2 ,T)6(s>2m 2 )- 



2s 



(27r) 3 " JV " '~ y ~ ~ "" '2E X 2E 2 
n 3 (T)((l - cos 19 )<Ty_nu (s > 2to 2 )) . 



(4) 



(5) 



where rij(T) is the thermal particle density in the rest 
frame of the QCD fluid, and where we have expressed 
the Mandelstam variable as s — 2E\E 2 {\ — cos 6*12), i.e. 
in terms of the angle B\ 2 between the hard parton and 
the thermal parton, and their energies in the fluid rest 
frame. Changing the integration variable cos #12 to s, it 
is straightforward to arrive at the following result 



ij—tkl 



1 f°° 

jg-ap / dE 2 f j (E 2 )Sl ij ^ kl (E 1 ,E 2 ,m 2 ), 



where 



(6) 



n iJ ^ k i(E 1 ,E 2 ,m 2 



AExE-2 



ds[saij^ki(s)}. (7) 



The analytical results for the relevant elements of f2y_>jy 
are listed in the Appendix. The remaining E 2 integral is 
easiest to do numerically. 

The obtained scattering rates Ti for quarks, antiquarks 
and gluons, along with the contributions ry_>.jy from dif- 
ferent processes, are shown in Fig. Q] as a function of tem- 
perature T and in Fig. [2] as a function of the hard parton 
energy E\ . 

To get a feeling of the general characteristics of the sys- 
tem, we consider also the following analytical estimates. 
In the process ij — > kl, the average of a quantity A in 
the plasma at a temperature T can be computed in the 
fluid rest frame as 



(A) 



1 



ij—¥kl 



f^-fJE 2 ,T)-A- 
(2tt) 3JA 2 ' ; ( 8 ) 

9(s > 2to 2 )(1 - cos#i2)er y ^fcz(s). 

At the high-energy limit, E 2 ,s >■ m 2 ,T 2 , the rele- 
vant pQCD cross sections all behave as a ~ mT 2 , and 

°gg^gg ~ I a gq^gq ~ (i) 2(T q,qj^g,qj ( see Appendix, take 
the limit E\ ^> to ~ E 2 , and s 3> \t\, \u\). In the Boltz- 
mann approximation we thus get 



(s)^8E 1 T, (E 2 )^3T, (cos0 12 ) 



(9) 
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• Total quark scattering rate 
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FIG. 1: (Color online) The scattering rates Vi of gluon (above) 
and quark (below) as a function of temperature T for a QCD 
plasma at rest. Flavour- and quark-antiquark -summed con- 
tributions from different processes and the analytical esti- 
mates for the total rates discussed in the text are also shown 
(solid lines). 



FIG. 2: (Color online) The scattering rates Ti of gluon (above) 
and quark (below) as a function of the parton energy E\ , for 
a QCD plasma at rest at T — 0.5 GeV. Flavour-summed con- 
tributions from different processes are shown. The analytical 
estimates for Ti (see text) are shown by the solid lines. In the 
upper panel, the analytical estimate for the gg — > gg process 
(dashed line) is the same as for the quark-antiquark summed 
99 99 process. In the lower panel, the dashed lines show the 
contributions from the qig — >• qig process and an individual 
qiqj — > qiqj process. 



from which we see that, in this limit, (s) = 2£i(£ , 2)(l — 
(cos #12)). The last result in Eq. © indicates that the 
collisions of the high-energy parton with the medium par- 
tons happen on the average in an angle 110°. 



Using the same approximation, we can easily obtain an- 
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alytical estimates also for the scattering rates, 



[r fi 



(g g + 2-3-g q -D- 
r 



9aT 3 



9 2 7TTO 2 



= 144- 



2 T 3 



(10) 



a? ff -i+2.3. s ,.(tn- 



9a 2 T 3 



2 7rm 2 



64 



a 2 T 3 



(11) 



where the factors 2 • 3 account for quarks and antiquarks 
times the number of quark flavours. These estimates are 
shown in Figs. [T] and Fig.[2J together with the actual rates 
applied in our simulation (once the boost to the local rest 
frame has been taken into account). We observe that the 
analytical rates describe the numerical ones (both the 
total rate and the ones from the main processes) quite 
well in the kinematic and temperature range of RHIC. 
In a static homogeneous QGP, the average number of 
collisions experienced by the high-energy parton of an 
average path length L, would, in the high-energy limit, 
be LT,. Taking L = 5 fm and using Eqs. (fTU1) - (fTT|) 
with a s = 0.3 and s m = 1, we arrive at the range 14— 
5 collisions for gluons and 6-2 collisions for quarks at 
T = 500 — 165 MeV. This gives us a simple order-of- 
magnitude estimate of the expected number of elastic 
scatterings for a hard parton. 



C. The Monte Carlo simulation 

To initiate the hard massless parton of a type i in each 
event, we sample the LO pQCD single-jet production 
spectrum (for more details, see [38|]) at yi — and at 

PTmin <PT < y/s/2, 



dcr 



AB^ri+X 



dyidy- 



E 



da 



AB^kl+X 



dp 1 T dy 1 dy 2 



[SkiS{y l - j/i ) 



S H S(yi -y 2 )} 



(12) 



where 



da 



dpfdy 1 dy 2 



AB^kl+X 



ah 



da ir 



->/.7 



Xlfa/A{%l,Q )X-Zh/ B {X2,Q 



(13) 

with do I dt containing the squared matrix elements listed 
in Eqs. (A-l)-(A-8). For the parton distribution func- 
tions (PDFs), we use the CTEQ6L1 set 3!) . For sim- 
plicity in this exploratory work we, however, neglect the 
nuclear effects to the PDFs [4CH42j | , since these are small 
in comparison with the ones arising from the final state 
interactions with the medium. The initial rapidity j/j is 
randomly generated in the range [y m in, Vmax] from a flat 
distribution. (This is consistent on the level of the lon- 
gitudinal boost-symmetry approximation made for the 



hydrodynamical part.) This fixes the hard-parton en- 
ergy E and polar angle 6 of its momentum vector. The 
azimuth angle 4> is evenly distributed between [0, 2ir]. 
In order to keep track of the hard partons interactions 
with the medium in the coordinate space, the initial spa- 
tial position of the hard parton has to be fixed as well. 
We start to follow the hard parton at the initial longitu- 
dinal proper time To of our hydrodynamical model. Since 
in the cm. frame all hard partons are produced in the 
Lorentz-contracted overlap region at z rs 0, the longitu- 
dinal position at later times (before the first collision at 
T > T o) is assumed to be determined by the longitudi- 
nal momentum only. The initial time and longitudinal 
coordinates for the hard parton are thus to = to cosh yi 
and zo = tq sinhj/i. The initial position in the transverse 
plane is sampled from the nuclear overlap function 



T AA (h)= J d 2 sT A {s + b/2)T A {s-b/2), 



(14) 



where b is the impact parameter. As we focus only on 
the central collisions, in the following b = 0. The nuclear 
thickness function T A (s) is defined as usual, 



Ta(s) 



dz Pa(s,z) 



(15) 



where for the nuclear density Pa(t) we use the standard 
Woods-Saxon distribution 



Pa{s,z) 



Pa 



expr^. 



1 



(16) 



with r = ^|s| 2 + z 2 , R A = 1.12AV3 - 0.86.4- 1 / 3 , d = 
0.54 fm and p = 0.17 fm~ 3 . 

The hard parton propagates through the plasma in small 
time steps At, during which we propagate the parton in 
position space. Since we consider only incoherent scat- 
terings here, the probability for not colliding in the time 
interval At is assumed to be given by the Poisson distri- 
bution 



P(No collisions in At) = e 



-r;At 



(17) 



where Ti is the total scattering rate ((TJ) for the hard par- 
ton of the type i. Hence the probability to collide at least 
once during the time At is l-e~ riAt sa T t +0((T t At) 2 ). 
For small enough At we can assume that there will be at 
most one collision. 

However, as explained above, we calculate the scattering 
rates (fTJ) in the local rest frame of the quark-gluon plasma 
fluid element and we must boost also the time step At 
to the same frame, 



At -> At' 



At 

5 

7 



where the Lorentz factor 7 = — ^=^= appears, u being 

the 3-dimensional flow velocity vector of the hydrody- 
namical plasma in the cm. frame of the colliding nuclei. 
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Should a scattering happen, the probability Pij^-ki for 
a given type of scattering process is determined by the 
partial scattering rates ((2]) as 



P(Process ij — > kl) 



r, 



(18) 



We produce the 4-momentum p2 of the scattering part- 
ner from the thermal medium according the energy and 
angular distributions defined by the differential scatter- 
ing rate dTi/d 3 p2 (see Eq. ((2J)) of the process. We then 
boost the 4-momenta of the hard parton pi and the ther- 
mal particle P2 to the cm. frame of the partonic collision, 
rotating the frame so that p\ is on the z-axis. It is now 
easy to determine the scattering angle 8* by sampling the 
distribution 



dt 



16irs 2 



using cos 8* — — + 1, and sampling the azimuthal angle 
ip* from a flat distribution. 

After determining the scattering angles 8* and ip*, we 
can construct the final state 4-momenta ps and p4 and 
Lorentz transform these back to the original cm. frame 
of the colliding nuclei. The final state parton with high- 
est energy is then chosen as the new hard parton to be 
propagated further, for which we repeat the procedure 
outlined above with the next timestep. Note that the 
type of the hard parton may change here. (We also note 
that in our code, we have the possibility to follow the 
recoil partons as well but at sufficiently high parton mo- 
menta (pr Si 5 GeV) , we do not see a significant effect 
on the observables studied here.) This procedure con- 
tinues until the surrounding matter has cooled down to 
the critical temperature Tq. At T = Tq the system is in 
this work in a mixed phase where energy density e keeps 
decreasing while the temperature stays constant. 
It is a priori not clear how the interaction with the 
medium in the mixed phase should be modelled properly 
in pQCD. However, we will demonstrate a posteriori (see 
Fig. O ahead) that within the assumption of an effective 
temperature, the effect of the mixed phase is small be- 
cause the medium density is small as compared to the ini- 
tial density. Thus, unless there is an anomalous enhance- 
ment of parton-medium scattering in the mixed phase by 
non-perturbative processes, a detailed treatment of the 
mixed phase is not important, and we can also ignore the 
interactions of the hard parton in the hadron gas. In the 
following, we treat the interactions in the mixed phase 
by using an effective temperature 



(M) 



where gQ — g g + ^2Nfg q = ^ is the quark-gluon plasma 
degrees of freedom and B is the bag constant. Our sim- 
ulation ends when T < Tq and the system enters the 
hadron gas phase, where we assume no significant inter- 
action between parton and medium. 



The outcome of the procedure described above is a 
medium-modified distribution of high-energy partons, 
dN dp T dy — ■ I n or der to compare with observables, we 
have to hadronize this distribution. In the current 
leading-particle case, this can be done by convoluting 
the obtained partonic distribution with the fragmenta- 
tion function Df^^(z, pb F ): 



dN 



AA^tt+X 



dP T dy 



dN AA^f+X 

dp T dy 



= X] / d Prdy 
f J 

/ dzDf^(z, n F )6(P T - zp T ), 
Jo 



(20) 



where z = — is the fraction of the final parton momen- 
turn pt available to the hadron with momentum Pt, and 
fip ~ Pt is the fragmentation scale. In this work we 
have used the Knichl-Kramcr-Pottcr (KKP) fragmenta- 
tion functions 1431 . 



III. RESULTS 

In the following, we choose a lower bound of the parton 
transverse momentum, pr min , as 5.0 GeV to account for 
the observation that at low Pt the main hadron produc- 
tion mechanism is not independent fragmentation. Par- 
tons which fall below E — 3T in our simulation are as- 
sumed to be thermalized and become part of the medium, 
and we do not follow these partons anymore. We average 
all observables across the rapidity window [—0.35, 0.35] 
which corresponds to the PHENIX acceptance. To make 
sure we account properly for all the possible partons 
falling into this final rapidity window, we choose the ini- 
tial y m in — —y ma x = — 1- We refer to all partons above 
both the PT min cut and within the final rapidity accep- 
tance as 'punch-through partons'. 

To verify the magnitude of the rapidity shifts experienced 
by the high-i? partons, we plot in Fig.[3]the initial rapid- 
ity distribution of the punch-through partons in our de- 
fault set-up with a s = 0.3 and s m = 1.0. We can see that 
about 13 % of the punch-through partons originate from 
outside of the final rapidity interval. On the one hand, 
this suggests that the eikonal approximation - where all 
partons would remain within the interval \y\ < 0.35 - is 
a fairly good first approximation. On the other hand, if 
one aims at a precision analysis of the single-hadron dis- 
tributions or, even more importantly, at correlations of 
two high-Pr hadrons, the non-eikonal dynamics should 
be obviously taken into account. An attempt to quanti- 
tatively understand the rapidity shift seen in Fig. [3] can 
be found in Appendix [Bl 

Also the leading-parton flavour-conversions are usually 
neglected in the eikonal energy-loss simulations. Inter- 
estingly, this effect is not symmetric between quarks and 
gluons: due the collision rates, there are more flavour 
conversions related to gluons than to quarks, and a quark 
from a g — > q leading-parton conversion is then more 



7 



In Fig.[4]we present the nuclear modification factor Raa, 
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FIG. 3: The initial rapidity distribution of the punch-through 
partons. The final interval \y\ < 0.35 is shown by the dashed 
lines. 



likely to become a punch-through parton than a gluon 
from ag->p conversion. Our simulation shows that for 
a simplified homogeneous static plasma cylinder of, say, 
radius L = 5 fm and T = 500 MeV, as much as 40 % of 
those original high-px gluons which will become punch- 
through partons, will turn into quarks. In the full simula- 
tion with 3-d hydrodynamical background medium, how- 
ever, only about 4% of all punch-through partons have 
undergone this effect. We thus conclude that in simu- 
lations for single-hadron distributions at RHIC, where a 
realistically evolving background medium is considered, 
the leading-parton flavour-conversions can be neglected. 
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dN\ A /dP T dy 
T A A(0)daPP/dP T dy' 



(21) 



for a s = 0.3, s m = 1.0, both on the hadronic level and on 
the partonic level, and compare it with the ir° data from 
PHENIX It is evident that for this 'default' choice 
of a s and s m elastic collisions are not enough to account 
for the data but that they should not be neglected, either. 
This result is in agreement with the calculation of purely 
elastic energy loss in the Arnold- Moore- Yaffe (AMY) for- 
malism averaged over hydrodynamics \'2 l! where a similar 
value of Raa ~ 0.6 is found for a s = 0.27. Note that a 
comparison with models not using a constrained fluid- 
dynamics model for the bulk evolution or the same a s is 
not meaningful. 

It can further be deduced from the figure that since the 
suppression as a function of Pt (or partonic pr) is flat, 
the convolution of the medium-modified parton spectrum 
with the fragmentation function does not cause a signif- 
icant effect in Raa- We will thus in the following focus 
on the partonic suppression pattern, with the implica- 
tion that the partonic suppression provides a reasonable 
estimate for the suppression also on the level of hadrons. 
As a technical point, we note that in order to collect 
enough statistics also at high-p^, the partonic results in 
Fig. d] consist of three dedicated runs, corresponding to 
p Tmm = 5 ; 10, 15 GeV. 
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FIG. 5: (Color online) The parton-level nuclear modification 
factors with our default set up a 3 = 0.3 and s m = 1.0, includ- 
ing only the scatterings in the QGP phase (squares) and in 
the QGP and mixed phase but not the in hadron gas (circles). 



FIG. 4: (Color online) Comparison of the tv° nuclear modi- 
fication factor both at hadron and parton level with strong 
coupling constant a s = 0.3. PHENIX data is from [44| . 



In Fig. [5] we demonstrate, as discussed above, that the 
mixed phase is in this set-up indeed causing only a small 
additional suppression and, consequently, that it is justi- 
fied to neglect the hadron gas effects in the present study. 
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The results shown in this and in the following figures are 
from runs with PT min = 5 GeV. 




i — ■ — i — ■ — i — ■ — i — ■ — i — ■ — r 
9 10 11 12 13 14 15 



p T [GeV] 



FIG. 6: (Color online) Parton level nuclear modification fac- 
tor Raa for the regulator-mass parameter value s m — 1.0 and 
strong coupling constant values a s — 0.15, 0.3, 0.6. 

In Fig. [6] we show the partonic nuclear modification for 
three different choices of constant a s . As expected, 
we observe a rather strong dependence on the value of 
a s . The behavior is quite straightforward; doubling the 
coupling approximately doubles the relative suppression. 
With a s = 0.6 we find a suppression due to elastic en- 
ergy loss already below the data. On the other hand, 
with a s = 0.15 the effect is rather weak. We thus focus 
on a s = 0.3 in the following when examining other model 
parameters. 

As is obvious the other significant model parameter in 
our simulation is s m , controlling the momentum cutoff 
scale m which regulates the scattering cross sections at 
smallest and largest angles. Like a s , it has a direct effect 
on the scattering cross sections of a hard parton with 
medium partons. However, it could be expected to in- 
fluence energy loss in a different way. While the total 
cross section can become very large if small-angle for- 
ward scattering processes are taken into account, such 
processes do not transfer a large amount of energy and 
are hence not efficient for energy loss. Thus, a smaller 
cutoff mass is expected to increase the total cross section 
without affecting large-angle scattering whereas a change 
in a s affects both forward and large-angle scattering in 
equal measure. 

A comparison of different values of the cutoff mass in 
Fig. [7] shows, however, that this parameter has an unex- 
pectedly large effect on the nuclear modification factor. 
The suppression increases with decreasing mass; this in- 
dicates that the cutoff mass is not yet in the region where 
forward scattering, which is inefficient for energy loss, 
takes place. 




10 11 12 13 14 15 



p T [GeV] 



FIG. 7: 

values s 
0.3. 



(Color online) Parton level Raa for mass parameter 
•n = 0.25, 0.5, 1.0 and strong coupling constant a s = 



In Fig. [5] we show the geometry underlying the sup- 
pression in terms of the initial spatial distribution of 
punch-through partons for different values of a s . This 
figure can be considered as a 'tomography' image, where 
we have rotated the position of each particle in trans- 
verse xy-plane in such a way that they all move to the 
same direction, to the direction of positive x-axis. With 
a s — 0.15 (upper panel), the quark-gluon plasma is quite 
transparent in transverse direction as far as elastic colli- 
sions are concerned, thus the initial distribution of punch- 
through partons is centered around the origin. Doubling 
the value of a s to 0.3 (middle panel) makes the distri- 
bution move off-center. Thus, the medium introduces a 
surface bias which means that hard partons from the cen- 
ter cannot punch through the matter anymore. Overall, 
the plasma becomes more opaque, as is expected from 
Fig. [SJ The opaqueness is emphasized even more in the 
bottom panel, which corresponds to a s — 0.6, which re- 
produces, or slightly overestimates the suppression mea- 
sured at RHIC. We have also checked the tomography 
images for different mass parameters (s m ) which repro- 
duce the same suppression pattern in Raa as when a s 
is varied (see Fig. [7]). Our conclusion here is that these 
do not significantly differ from the ones we have shown 
in Fig. 1 



IV. DISCUSSION 

The results above indicate on the one hand that for our 
default choice of parameters a s — 0.3 and s m — 1, elas- 
tic energy loss is not sufficient to account for the ob- 
served value of Raa- On the other hand, we see that 
at least in the framework which is based on individual 
parton-parton scatterings instead of, e.g., static scatter- 
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account for the physics of a dynamically running a s and 
thermal field-theoretical processes generating a screening 
mass. Thus for a reasonable variation of these parame- 
ters there are pairs (a s , s m ) for which the Raa data can 
be described by elastic processes alone. It is thus evi- 
dent from these findings that a comparison based on the 
magnitude of the single hadron suppression alone cannot 
resolve the question how significant the contribution of 
2 — !• 2 processes to energy loss actually is. To improve 
upon this problem, the next step is to include 2^3 
processes and coherence effects into the simulation. 

However, a comparison of the tomographic pictures in 
this paper (Fig. [3 with the results of a radiative energy 
loss calculation [12| indicates that the way hard partons 
probe the medium geometry through energy loss is quite 
different for radiative and elastic energy loss. Thus, as 
outlined in a schematical way in [25] it can be expected 
that pathlength-dependent observables, like the suppres- 
sion as a function of the angle of the observed hadron 
with the reaction plane, or the suppression of back-to- 
back correlations, will show pronounced differences be- 
tween elastic and radiative energy loss. The MC code 
introduced here provides an excellent tool to explore this 
idea, and this will be addressed in a future publication. 
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FIG. 8: (Color online) Initial production points of the punch- 
through partons on the transverse plane for a 3 = 0.15 (above) 
a s = 0.3 (middle )and a s — 0.6 (below). All particles are 
moving to the positive ^-direction. 



ing centers, the elastic energy losses give a non-negligible 
contribution to Raa- This has also been observed by a 
comparable calculation within the AMY formalism [24[ in 
which a combination of radiative and elastic energy loss 
was needed to account for the data. However, our re- 
sults show an unfortunately strong dependence on cutoff 
parameters which are in essence supposed to effectively 



Appendix A: Analytical expressions for scattering 
rates 



The spin- and color summed/averaged squared parton 
scattering amj 
(see e.g. [371 ]) 



scattering amplitudes IMj^-^^ for gluons and quarks are 



ij — > kl 



\M\ 2 



ij—¥kl 



99 -> 99 

99 -> qq 

gq -> gq 
gq -> gq 

q-q, -> '/.'/.. 
qiq 2 -> <Mj 
Mj -> <Mj 

qtqi -> 
<Mi -> q%qi 

q,q. -> '/..'/.. 

qm qm 

qq -> 99 
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(A-3) 
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Inserting the expression §4$ for the cross section into the 
definition of the ^-function ([7]) we have 



^ij^kl(Ei,E 2 ,m ) 



1 + Ski 167T J 2 



dt\M l2 



s-\-m 2 



ij—^-kl 



Several amplitudes contain identical terms. We list the 
integrals for each of these terms separately, as it is then 
easy to combine the appropriate terms to form the final 
expressions for Uij^ki (E\ , E 2 , m 2 ) . 



X(s,t,u) 



&E1E2 
2m 2 



ds 



i f~™ 2 dtX(s,t,u) 

s J —s-\-m 2 V J J / 



4EiE 2 - 2m 2 log 2 ^ 2 - - 2m 2 
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ut 



su st 
t 2 +u 2 
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(A-10) 
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(A-12) 
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Appendix B: Rapidity shift estimates 

Here we try to understand quantitatively the magnitude 
of the rapidity shifts of punch-through partons shown in 
Fig. |3l At the high-energy limit E\ 3> T,m, a straight- 
forward calculation in the static plasma approximation 
(see Sec. IIIBj) gives the average scattering angle in the 
CMS frame as 

(cosfl*) = 1 + (2t/a) « 1 - log(^). (B-l) 

We set the initial high-i? parton at zero rapidity and 
its momentum pi to point into the direction of the x- 
axis, while the beam direction is, as usual, the z-axis. 
Having the maximum rapidity shift in mind, we take the 
medium-parton's momentum p 2 to reside in the (x, z)- 
plane, its energy to be E 2 — 3T and the collision (polar) 
angle relative to pi such that cos#i2 = —1/3, which 
correspond to the averages in Eq. |9] Thus 

Pi = e x i?i, p 2 = e x E 2 cos Q 12 + e z E 2 sin 9 12 . (B-2) 

In a reverse boost from the CMS frame of the parton- 
parton scattering, the energy and momentum of the scat- 
tered high-E parton become 

E 3 = 7 (£*+vcms -Ps*) (B-3) 
P3 = P3*+7 v Cms( — — rP3* • V C MS +E*),(B-A) 



where the velocity of the CMS frame is vqms — (Pi + 
Pi)/(Ei + E 2 ), 7 is the Lorentz gamma factor and the 
asterisks denote the CMS-frame quantities. The momen- 
tum of the scattered high-i? parton in the CMS-frame 
can be computed to be 

Pa* = e x [-p^sin^sin^ +pJ !B cos0£] 

+ e y E* sin 9% cos (f>% (B-5) 
+ e z [p* x sinf|sin^3 + pl z cos 9%}, 

where #3 is the scattering (polar) angle relative to pi* 
and $3 is the azimuthal angle around pi * , measured from 
the y-axis, and where the CMS-frame quantities are ob- 
tained in a boost to vcms as 

E* = 7(^1 -vcms -Pi) (B-6) 

7 

Pi* = Pi+7 v Cms( — ; — 7P1 ' v cms — Ei). (B-7) 
7+1 

The point here is that although in the CMS frame (where 
P2* = — Pi*) the scattering cross section does not de- 
pend on the angle <j>*, a 0* dependence is generated in 
the boost back to the plasma rest frame. The rapidity 
shift Sy — 2/3 is now easy to obtain numerically from 
the components E3 and p 3z and we find that to a good 
approximation 

8y & Asmfe. (B-8) 

For a hard parton of, say, Ei = 10 GeV in a T = 500 MeV 
plasma, we obtain A w 0.26. The maximum shifts 
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FIG. 9: (Color online) Rapidity spread of punch-through par- 
tons which originate from a 10 GeV, y = parton at the 
spatial origin. 



Sy ss ±A occur at 4>* = ir/2 and 37r/2, i.e. when the 
scattering remains in the (x, z)-plane. Similarly, we find 
that E 3 /GcV fa 8.7+ 1.4 sin 0* and 6 3 « 15° + 1.1° sin^. 
To estimate the maximum possible total rapidity shift 
|AY| induced by n scatterings, we may assume that at 
the limit of small scattering angles (small rapidity shifts) 
the maximum-shift distribution Sy rj A sin (jf holds in 
each scattering. The rapidity shifts are additive and we 
can write the probability distribution for AY as follows: 

„ n 
P(AY)= / d6y 1 ...d6y n P(6y 1 )...P(6y n )5(AY-J2Syi), 

(B-9) 

where the probability distribution for a rapidity shift in 
one collision, 



For an order of magnitude estimate, we take n ~ 6 (from 
our simulation, since for this static high-T plasma test 
case the leading-parton flavour-coversions are in fact a 
non-negligible effect), and arrive at the average maximum 
total rapidity shift |AY| = ^((AY) 2 ) = A^/^/2 : 



0.45. 

Figure [5] shows our simulation test results for the ra- 
pidity spread of punch-through partons which originate 
from 10 GeV partons sent off from the origin into the 
x-direction. The width of the distribution in the static 
500 MeV plasma case (squares) is seen to be about 0.3, 
which nicely confirms the order of magnitude of our up- 
per limit estimate |AY|. The smaller width in the simula- 
tion follows from the fact that the simulation accounts for 
scatterings (corresponding to a continuum of azimuthal 
incoming angles <f>i2 around pi), where the rapidity shift 
effects are smaller than in our estimate above. Changing 
the background into a hydrodynamically evolving more 
realistic medium (circles) reduces the width of the ra- 
pidity shift by about a factor 3. This is caused by the 
fact that collision rates, scattering angles #3 and rapidity 
kicks all decrease with decreasing temperature. Finally, 
in a full simulation, where the high-E partons can be 
produced anywhere in the transverse plane, and where 
also the surface bias effect is present (see Fig. [8]) , the ra- 
pidity shifts become still somewhat smaller. From these 
considerations, and realizing that the rapidity spread dis- 
tributions should be applied to each final-rapidity bin of 
Fig. [3j we conclude that the degree of non-eikonality (y- 
shifts) is understood, and also that if we wish to fully 
account for punch-through partons in the finite rapidity 
interval \y\ < 0.35, we have to consider an initial window 
of j^j < 1. " 
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P(Sy) = 



(B-10) 



where \5y\ < A, is obtained from the corresponding flat 
azimuthal angle distribution P(4>*) = l/(27r) by a change 
of variables according to Eq. IB-81 Now it is easy to obtain 

((AY) 2 ) = n((6y) 2 } = n^-. (B-ll) 
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